
      module mo_strato_rates
!=======================================================================
! ROUTINE
!   ratecon_sfstrat.f
!
!  Date...
!  15 August 2002
!  11 April  2008
!
!  Programmed by...
!   Douglas E. Kinnison
!
! DESCRIPTION
!
! Derivation of the rate constant for reactions on
!   sulfate, NAT, and ICE aerosols.
!
!
! Sulfate Aerosol Reactions               Rxn#   Gamma
!   N2O5   + H2O(l)     =>  2HNO3         (1)    f(wt%)
!   ClONO2 + H2O(l)     =>  HOCl + HNO3   (2)    f(T,P,HCl,H2O,r)
!   BrONO2 + H2O(l)     =>  HOBr + HNO3   (3)    f(T,P,H2O,r)
!   ClONO2 + HCl(l)     =>  Cl2  + HNO3   (4)    f(T,P,HCl,H2O,r)
!   HOCl   + HCl(l)     =>  Cl2  + H2O    (5)    f(T,P,HCl,HOCl,H2O,r)
!   HOBr   + HCl(l)     =>  BrCl + H2O    (6)    f(T,P,HCl,HOBr,H2O,r)
!
! Nitric Acid Di-hydrate Reactions        Rxn#    Gamma   Reference
!   N2O5   + H2O(s)     =>  2HNO3         (7)     4e-4   JPL06-2
!   ClONO2 + H2O(s)     =>  HOCl + HNO3   (8)     4e-3   JPL06-2
!   ClONO2 + HCl(s)     =>  Cl2  + HNO3   (9)     0.2    JPL06-2
!   HOCl   + HCl(s)     =>  Cl2  + H2O    (10)    0.1    JPL06-2
!   BrONO2 + H2O(s)     =>  HOBr + HNO3   (11)    0.3    David Hanson PC
!
! ICE Aersol Reactions                    Rxn#    Gamma
!   N2O5   + H2O(s)     =>  2HNO3         (12)     0.02   JPL06-2
!   ClONO2 + H2O(s)     =>  HOCl + HNO3   (13)     0.3    JPL06-2
!   BrONO2 + H2O(s)     =>  HOBr + HNO3   (14)     0.3    JPL06-2
!   ClONO2 + HCl(s)     =>  Cl2  + HNO3   (15)     0.3    JPL06-2
!   HOCl   + HCl(s)     =>  Cl2  + H2O    (16)     0.2    JPL06-2
!   HOBr   + HCl(s)     =>  BrCl + H2O    (17)     0.3    JPL06-2
!
! NOTE: The rate constants derived from species reacting with H2O are
!       first order (i.e., sec-1 units) - an example is N2O5 + H2O = 2HNO3.
!       Other reactions, e.g., ClONO2 + HCl have rate constants that
!       are second order (i.e., cm+3 molecules-1 sec-1 units). In all
!       of these types of reactions the derived first order rate constant
!       {0.25*(mean Velocity)*SAD*gamma} is divided by the HCl abundance
!       to derive the correct second order units.
!
! NOTE: Liquid Sulfate Aerosols...
!       See coding for references on how the Sulfate Aerosols were handled.
!       Data was used that was more recent than JPL00.
!
!
! INPUT:
!  ad      .    .... air density, molec. cm-3
!  pmid        ..... pressures, hPa
!  temp        ..... temperatures, K
!  rad_sulfate ..... Surface area density, cm2 cm-3
!  sad_sulfate ..... Surface area density, cm2 cm-3
!  sad_nat     ..... Surface area density, cm2 cm-3
!  sad_ice     ..... Surface area density, cm2 cm-3
!  brono2mv    ..... BrONO2 Volume Mixing Ratio
!  clono2mvr   ..... ClONO2 Volume Mixing Ratio
!  h2omvr      ..... H2O Volume Mixing Ratio
!  hclmvr      ..... HCl Volume Mixing Ratio
!  hobrmvr     ..... HOBr Volume Mixing Ratio
!  hoclmvr     ..... HOCl Volume Mixing Ratio
!  n2o5mvr     ..... N2O5 Volume Mixing Ratio
!
! OUTPUT:
!
!  rxt         ..... Rate constant (s-1 and cm3 sec-1 molec-1)
!=======================================================================

      private
      public :: ratecon_sfstrat, init_strato_rates, has_strato_chem

      integer :: id_brono2, id_clono2, id_hcl, id_hocl, &
           id_hobr, id_n2o5
      integer :: rid_het1,  rid_het2,  rid_het3,  rid_het4,  rid_het5, &
           rid_het6,  rid_het7,  rid_het8,  rid_het9,  rid_het10, &
           rid_het11, rid_het12, rid_het13, rid_het14, rid_het15, &
           rid_het16, rid_het17

      logical :: has_strato_chem 

      contains

        subroutine init_strato_rates

          use mo_chem_utls, only : get_rxt_ndx, get_spc_ndx
          use mo_aero_settling, only: strat_aer_settl_init
          implicit none

          integer :: ids(23)

          rid_het1  = get_rxt_ndx( 'het1' )
          rid_het2  = get_rxt_ndx( 'het2' )
          rid_het3  = get_rxt_ndx( 'het3' )
          rid_het4  = get_rxt_ndx( 'het4' )
          rid_het5  = get_rxt_ndx( 'het5' )
          rid_het6  = get_rxt_ndx( 'het6' )
          rid_het7  = get_rxt_ndx( 'het7' )
          rid_het8  = get_rxt_ndx( 'het8' )
          rid_het9  = get_rxt_ndx( 'het9' )
          rid_het10 = get_rxt_ndx( 'het10' )
          rid_het11 = get_rxt_ndx( 'het11' )
          rid_het12 = get_rxt_ndx( 'het12' )
          rid_het13 = get_rxt_ndx( 'het13' )
          rid_het14 = get_rxt_ndx( 'het14' )
          rid_het15 = get_rxt_ndx( 'het15' )
          rid_het16 = get_rxt_ndx( 'het16' )
          rid_het17 = get_rxt_ndx( 'het17' )

          id_brono2 = get_spc_ndx( 'BRONO2' )
          id_clono2 = get_spc_ndx( 'CLONO2' )
          id_hcl    = get_spc_ndx( 'HCL' )
          id_hocl   = get_spc_ndx( 'HOCL' )
          id_hobr   = get_spc_ndx( 'HOBR' )
          id_n2o5   = get_spc_ndx( 'N2O5' )

          ids(:) = (/ rid_het1, rid_het2, rid_het3, rid_het4, rid_het5, rid_het6, rid_het7, rid_het8, &
               rid_het9, rid_het10, rid_het11, rid_het12, rid_het13, rid_het14, rid_het15, &
               rid_het16, rid_het17, id_brono2, id_clono2, id_hcl, id_hocl, id_hobr, id_n2o5 /)

          has_strato_chem = all( ids(:) > 0 )

          if (.not. has_strato_chem) return

          call strat_aer_settl_init

        endsubroutine init_strato_rates

      subroutine ratecon_sfstrat( ad, pmid, temp, rad_sulfate, sad_sulfate, &
                                  sad_nat, sad_ice, h2ovmr, vmr, rxt, ncol, &
                                  aoc, aso4, wtper, wtorg, mh2so4, &
								  gprob_n2o5, &
								  gprob_cnt_hcl, &
								  gprob_cnt_h2o, &
								  gprob_bnt_h2o, &
								  gprob_hocl_hcl, &
								  gprob_hobr_hcl, &
								  Hhcl,&
								  Mhcl,&
								  rlats)

      use shr_kind_mod, only : r8 => shr_kind_r8
      use chem_mods,    only : adv_mass, rxntot, gas_pcnst
      use ppgrid,       only : pcols, pver
      use mo_sad,       only : sad_top
      use cam_logfile,  only : iulog

      implicit none

!-----------------------------------------------------------------------
!	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in) :: ncol                                 ! columns in chunk
      real(r8), dimension(ncol,pver,gas_pcnst), intent(in) :: &   ! species concentrations (mol/mol)
        vmr
      real(r8), dimension(ncol,pver), intent(in) :: &
        ad, &                                                     ! Air Density (molec. cm-3)
        rad_sulfate, &                                            ! Radius of Sulfate Aerosol (cm)
        sad_ice, &                                                ! ICE Surface Area Density (cm-1)
        sad_nat, &                                                ! NAT Surface Area Density (cm-1)
        sad_sulfate, &                                            ! Sulfate Surface Area Density (cm-1)
        h2ovmr                                                    ! water vapor volume mixing ratio( gas phase )
      real(r8), dimension(pcols,pver), intent(in) :: &
        pmid, &                                                   ! pressure (Pa)
        temp                                                      ! temperature (K)

      real(r8), intent(out) :: &
        rxt(ncol,pver,rxntot)                                     ! rate constants
        
      ! Added for CARMA organics correction for weight percent
      real(r8), intent(in), optional :: &
        aoc(ncol,pver), &
        aso4(ncol,pver) 

	  real(r8), dimension(ncol,pver), intent(out) :: &            ! diagnostics
        gprob_n2o5, &
        gprob_cnt_hcl, &
        gprob_cnt_h2o, &
        gprob_bnt_h2o, &
        gprob_hocl_hcl, &
        gprob_hobr_hcl, &
        wtper, wtorg, mh2so4, Hhcl, Mhcl
        
       real(r8), intent(in) :: rlats(ncol) ! latitude in radiance
        
!-----------------------------------------------------------------------
!  	... local variables
!-----------------------------------------------------------------------
	real(r8), parameter :: small_conc = 1.e-30_r8
	real(r8), parameter :: av_const   = 2.117265e4_r8  ! (8*8.31448*1000 / PI)
	real(r8), parameter :: pa2mb      = 1.e-2_r8       ! Pa to mb
	real(r8), parameter :: m2cm       = 100._r8        ! meters to cms

      integer :: &
        i, &                      ! altitude loop index
        k, &                      ! level loop index
        m                         ! species index
        
      real(r8) :: wt_org, wt_water

!-----------------------------------------------------------------------
!   	... variables for gamma calculations
!-----------------------------------------------------------------------
      real(r8) :: &
        brono2vmr, &                            ! BrONO2 Volume Mixing Ratio
        clono2vmr, &                            ! ClONO2 Volume Mixing Ratio
        hclvmr, &                               ! HCl Volume Mixing Ratio
        hcldeni, &                              ! inverse of HCl density
        cntdeni, &                              ! inverse of ClONO2 density
        hocldeni, &                             ! inverse of HOCl density
        hobrdeni, &                             ! inverse of HOBr density
        hoclvmr, &                              ! HOCl Volume Mixing Ratio
        hobrvmr, &                              ! HOBr Volume Mixing Ratio
        n2o5vmr                                 ! N2O5 Volume Mixing Ratio

        real(r8) :: &
        av_n2o5, &                              ! N2O5 Mean Velocity (cm s-1)
        av_clono2, &                            ! ClONO2 Mean Velocity (cm s-1)
        av_brono2, &                            ! BrONO2Mean Velocity (cm s-1)
        av_hocl, &                              ! HOCl Mean Velocity (cm s-1)
        av_hobr                                 ! HOBr Mean Velocity (cm s-1)

      real(r8) :: &
        pzero_h2o, &                            ! H2O sat vapor press (mbar)
        e0, e1, e2, e3, &                       ! coefficients for H2O sat vapor press.
        aw, &                                   ! Water activity
        m_h2so4, &                              ! H2SO4 molality (mol/kg)
        wt, &                                   ! wt % H2SO4
        y1, y2, &                               ! used in H2SO4 molality
        &  a1, b1, c1, d1, &                    ! used in H2SO4 molality
        a2, b2, c2, d2                          ! used in H2SO4 molality

        real(r8) :: &
        z1, z2, z3, &                           ! used in H2SO4 soln density
        den_h2so4, &                            ! H2SO4 soln density, g/cm3
        mol_h2so4, &                            ! Molality of H2SO4, mol / kg
        molar_h2so4, &                          ! Molarity of H2SO4, mol / l
        x_h2so4, &                              ! H2SO4 mole fraction
        x_hcl, &                                ! HCl mole fraction
        ratio_hcl, &                            ! HCl mole ratio 
        den_hex, &                              ! density of hexanoic acid (g/L)
        aconst, tzero, &                        ! used in viscosity of H2SO4
        vis_h2so4, &                            ! H2SO4 viscosity
        ah, &                                   ! Acid activity, molarity units
        term1,term2,term3,term4, &              ! used in ah
        term5,term6,term7,term0, &
        T_limit, &                              ! temporary variable for temp (185-260K range)
        T_limiti, &                             ! 1./T_limit
        T_limitsq, &                            ! sqrt( T_limit )
        rad_sulf, &                             ! temporary variable for sulfate radius (cm)
        sadsulf, &                              ! temporary variable for sulfate radius (cm)
        sadice, &                               ! temporary variable for sulfate radius (cm)
        sadnat                                  ! temporary variable for sulfate radius (cm)

      real(r8) :: &
        C_cnt, S_cnt, &                         ! used in H_cnt
        H_cnt, &                                ! Henry's law coeff. for ClONO2
        H_hcl, &                                ! Henry's law coeff. for HCl
        D_cnt, &
        k_hydr, &
        k_h2o, &
        k_h, &
        k_hcl, &
        rdl_cnt, &
        f_cnt, &
        M_hcl, &
        atmos

      real(r8) :: &
        Gamma_b_h2o, &
        Gamma_cnt_rxn, &
        Gamma_b_hcl, &
        Gamma_s, &
        Fhcl, &
        Gamma_s_prime, &
        Gamma_b_hcl_prime, &
        Gamma_b, &
!        gprob_n2o5, &
        gprob_rxn, &
        gprob_tot, &
        gprob_cnt!, &
!        gprob_cnt_hcl, &
!        gprob_cnt_h2o

        real(r8) :: &
        D_hocl, &
        k_hocl_hcl, &
        C_hocl, &
        S_hocl, &
        H_hocl, &
        Gamma_hocl_rxn, &
        rdl_hocl, &
        f_hocl!, &
!        gprob_hocl_hcl

        real(r8) :: &
        h1, h2, h3, &
        alpha!, &
!        gprob_bnt_h2o

      real(r8) :: &
        C_hobr, &
        D_hobr, &
        aa, bb, cc, dd, &
        k_hobr_hcl, &
        k_dl, &
        k_wasch, &
        H_hobr, &
        rdl_hobr, &
        Gamma_hobr_rxn, &
        f_hobr!, &
!        gprob_hobr_hcl

      real(r8) :: &
        pmb,&					! Pressure, mbar (hPa)
        pH2O_atm,&				! Partial press. H2O (atm)
        pH2O_hPa,&				! Partial press. H2O (hPa)
        pHCl_atm,&				! Partial press. HCl (atm)
        pCNT_atm                                ! Partial press. ClONO2 (atm)

!-----------------------------------------------------------------------
!     	... Used in pzero h2o calculation
!-----------------------------------------------------------------------
      real(r8), parameter :: wt_e0 = 18.452406985_r8
      real(r8), parameter :: wt_e1 = 3505.1578807_r8
      real(r8), parameter :: wt_e2 = 330918.55082_r8
      real(r8), parameter :: wt_e3 = 12725068.262_r8

      real(r8) :: &
        wrk, tmp

      real(r8), parameter :: small = 1.e-16_r8

      if (.not. has_strato_chem) return

!-----------------------------------------------------------------------
!     	... intialize rate constants
!-----------------------------------------------------------------------
      do k = 1,pver
         rxt(:,k,rid_het1) = 0._r8
         rxt(:,k,rid_het2) = 0._r8
         rxt(:,k,rid_het3) = 0._r8
         rxt(:,k,rid_het4) = 0._r8
         rxt(:,k,rid_het5) = 0._r8
         rxt(:,k,rid_het6) = 0._r8
         rxt(:,k,rid_het7) = 0._r8
         rxt(:,k,rid_het8) = 0._r8
         rxt(:,k,rid_het9) = 0._r8
         rxt(:,k,rid_het10) = 0._r8
         rxt(:,k,rid_het11) = 0._r8
         rxt(:,k,rid_het12) = 0._r8
         rxt(:,k,rid_het13) = 0._r8
         rxt(:,k,rid_het14) = 0._r8
         rxt(:,k,rid_het15) = 0._r8
         rxt(:,k,rid_het16) = 0._r8
         rxt(:,k,rid_het17) = 0._r8
      end do

	  do k = 1, pver
		 gprob_n2o5(:,k)      = 0._r8
		 gprob_cnt_h2o(:,k)   = 0._r8
		 gprob_cnt_hcl(:,k)   = 0._r8
		 gprob_bnt_h2o(:,k)   = 0._r8
		 gprob_hocl_hcl(:,k)  = 0._r8
		 gprob_hobr_hcl(:,k)  = 0._r8
		 wtper(:,k)           = 0._r8
         wtorg(:,k)           = 0._r8
         mh2so4(:,k)          = 0._r8
         Hhcl(:,k)            = 0._r8
         Mhcl(:,k)            = 0._r8
                  
	  end do

!-----------------------------------------------------------------------
!     	... set rate constants
!-----------------------------------------------------------------------
Level_loop : &
      do k = sad_top+1,pver
column_loop : &
         do i = 1,ncol
!-----------------------------------------------------------------------
!	... set species, pmb, and atmos
!-----------------------------------------------------------------------
	    brono2vmr = vmr(i,k,id_brono2)
	    clono2vmr = vmr(i,k,id_clono2)
	    hclvmr    = vmr(i,k,id_hcl)
	    hoclvmr      = vmr(i,k,id_hocl)
	    hobrvmr      = vmr(i,k,id_hobr)
	    if( hclvmr > 0._r8 ) then
	       hcldeni  = 1._r8/(hclvmr*ad(i,k))
	    end if
	    if( clono2vmr > 0._r8 ) then
	       cntdeni  = 1._r8/(clono2vmr*ad(i,k))
	    end if
	    if( hoclvmr > 0._r8 ) then
	       hocldeni  = 1._r8/(hoclvmr*ad(i,k))
	    end if
	    if( hobrvmr > 0._r8 ) then
	       hobrdeni  = 1._r8/(hobrvmr*ad(i,k))
	    end if
	    n2o5vmr      = vmr(i,k,id_n2o5)
            sadsulf      = sad_sulfate(i,k)
            sadnat       = sad_nat(i,k)
	    sadice       = sad_ice(i,k)
            pmb          = pa2mb*pmid(i,k)
            atmos        = pmb/1013.25_r8

!-----------------------------------------------------------------------
!  	... setup for stratospheric aerosols
!           data range set: 185K - 240K;    GRL, 24, 1931, 1997
!-----------------------------------------------------------------------
            T_limit   = max( temp(i,k),185._r8 )
            T_limit   = min( T_limit,240._r8 )
            T_limiti  = 1._r8/T_limit
            T_limitsq = sqrt( T_limit )

!-----------------------------------------------------------------------
!     .... Average velocity (8RT*1000/(PI*MW))**1/2 * 100.(units cm s-1)
!     .... or (av_const*T/M2)**1/2
!-----------------------------------------------------------------------
	    wrk       = av_const*T_limit
            av_n2o5   = sqrt( wrk/adv_mass(id_n2o5) )*m2cm
            av_clono2 = sqrt( wrk/adv_mass(id_clono2) )*m2cm
            av_brono2 = sqrt( wrk/adv_mass(id_brono2) )*m2cm
            av_hocl   = sqrt( wrk/adv_mass(id_hocl) )*m2cm
            av_hobr   = sqrt( wrk/adv_mass(id_hobr) )*m2cm
has_sadsulf : &
            if( sadsulf > 0._r8 ) then
!-----------------------------------------------------------------------
!     .... Partial Pressure of H2O, ClONO2, and HCl in atmospheres
!-----------------------------------------------------------------------
               if( hclvmr > 0._r8 ) then
                  pHCl_atm  = hclvmr*atmos
               else
                  pHCl_atm  = 0._r8
               end if

               if( clono2vmr > 0._r8 ) then
                  pCNT_atm  = clono2vmr*atmos
               else
                  pCNT_atm  = 0._r8
               end if

               if( h2ovmr(i,k) > 0._r8 ) then
                  pH2O_atm  = h2ovmr(i,k)*atmos
               else
                  pH2O_atm  = 0._r8
               end if

!-----------------------------------------------------------------------
!     .... Partial Pressure of H2O in hPa
!-----------------------------------------------------------------------
               pH2O_hpa = h2ovmr(i,k)*pmb
!-----------------------------------------------------------------------
!     .... Calculate the h2so4 Wt% and Activity of H2O - 
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     ... Saturation Water Vapor Pressure (mbar)
!-----------------------------------------------------------------------
               pzero_h2o = exp( wt_e0 - T_limiti*(wt_e1 + T_limiti*(wt_e2 - T_limiti*wt_e3)) )
!-----------------------------------------------------------------------
!     ... H2O activity
!     ... if the activity of H2O goes above 1.0, wt% can go negative
!-----------------------------------------------------------------------
               aw = ph2o_hpa / pzero_h2o
               aw = min( aw,1._r8 )
               aw = max( aw,.001_r8 )
!-----------------------------------------------------------------------
!     ... h2so4 Molality (mol/kg)
!-----------------------------------------------------------------------
               if( aw <= .05_r8 ) then
                  a1 = 12.37208932_r8
                  b1 = -0.16125516114_r8
                  c1 = -30.490657554_r8
                  d1 = -2.1133114241_r8
                  a2 = 13.455394705_r8
                  b2 = -0.1921312255_r8
                  c2 = -34.285174607_r8
                  d2 = -1.7620073078_r8
               else if( aw > .05_r8 .and. aw < .85_r8 ) then
                  a1 = 11.820654354_r8
                  b1 = -0.20786404244_r8
                  c1 = -4.807306373_r8
                  d1 = -5.1727540348_r8
                  a2 = 12.891938068_r8
                  b2 = -0.23233847708_r8
                  c2 = -6.4261237757_r8
                  d2 = -4.9005471319_r8
               else
                  a1 = -180.06541028_r8
                  b1 = -0.38601102592_r8
                  c1 = -93.317846778_r8
                  d1 = 273.88132245_r8
                  a2 = -176.95814097_r8
                  b2 = -0.36257048154_r8
                  c2 = -90.469744201_r8
                  d2 = 267.45509988_r8
               end if
!-----------------------------------------------------------------------
!     ... h2so4 mole fraction
!-----------------------------------------------------------------------
               y1       = a1*(aw**b1) + c1*aw + d1
               y2       = a2*(aw**b2) + c2*aw + d2
               m_h2so4  = y1 + ((T_limit - 190._r8)*(y2 - y1)) / 70._r8
               mh2so4(i,k) = m_h2so4
!-----------------------------------------------------------------------
!     ... h2so4 Weight Percent
!-----------------------------------------------------------------------
               wt = 9800._r8*m_h2so4 / (98._r8*m_h2so4  + 1000._r8)

              ! work with Susan Solomon and Doug Kinnison
              ! pyrocb organics dilutes aerosol solution, makes h2so4 weight
              ! percent smaller
              ! Feb.2 2022
              ! Pengfei Yu
              ! not all organics can be accounted,
              ! apply a simple mass ratio here
              ! Feb.19 2022
    
              wtorg(i,k) = wt ! sulfate weight percent without considering organics 
    
              if (present(aoc) .and. present(aso4) .and. aso4(i,k) .gt. 0._r8) then
               wt = 1._r8 / (1._r8 / wt + aoc(i,k) / aso4(i,k) / 100._r8)
               wt = max(wt, 0.1_r8)
              end if
              
              wtper(i,k) = wt ! sulfate weight percent considering organics
    
              ! Kane Stone April 12 2022
              ! Updating to include wt_org and wt_water for h2so4 mole fraction calculation  
              wt_org = wt * aoc(i,k) / aso4(i,k)
              wt_water = 100._r8 - wt - wt_org

!-----------------------------------------------------------------------
!     .... Parameters for h2so4 Solution, JPL-00
!-----------------------------------------------------------------------
!     ... h2so4 Solution Density (g/cm3)
!-----------------------------------------------------------------------
	       wrk = T_limit*T_limit
               z1 =  .12364_r8  - 5.6e-7_r8*wrk
               z2 = -.02954_r8  + 1.814e-7_r8*wrk
               z3 =  2.343e-3_r8 - T_limit*1.487e-6_r8 - 1.324e-8_r8*wrk
!-----------------------------------------------------------------------
!     ... where mol_h2so4 is molality in mol/kg
!-----------------------------------------------------------------------
               den_h2so4 = 1._r8 + m_h2so4*(z1 + z2*sqrt(m_h2so4) + z3*m_h2so4)
!-----------------------------------------------------------------------
!     ... h2so4 Molarity, mol / l
!-----------------------------------------------------------------------
               molar_h2so4 = den_h2so4*wt/9.8_r8
!-----------------------------------------------------------------------
!     ... h2so4 Mole fraction
!-----------------------------------------------------------------------
               ! Kane Stone April 12 2022
               ! Updating x_h2so4 calculation to account for solution being part water (molar mass 18) and part organic acid (molar
               ! mass 116 for hexanoic) 

               x_h2so4   = wt / (wt + (wt_water * 98._r8/18._r8) + (wt_org * 98._r8 / 116._r8))
               
               term1     = .094_r8 - x_h2so4*(.61_r8 - 1.2_r8*x_h2so4)
               term2     = (8515._r8 - 10718._r8*(x_h2so4**.7_r8))*T_limiti
               H_hcl     = term1 * exp( -8.68_r8 + term2 )
               M_hcl     = H_hcl*pHCl_atm
               
               Hhcl(i,k) = H_hcl
               Mhcl(i,k) = M_hcl

!-----------------------------------------------------------------------
!     ... h2so4 solution viscosity
!-----------------------------------------------------------------------
               aconst    = 169.5_r8 + wt*(5.18_r8 - wt*(.0825_r8 - 3.27e-3_r8*wt))
               tzero     = 144.11_r8 + wt*(.166_r8 - wt*(.015_r8 - 2.18e-4_r8*wt))
               vis_h2so4 = aconst/(T_limit**1.43_r8) * exp( 448._r8/(T_limit - tzero) )

!-----------------------------------------------------------------------
!     ... Acid activity in molarity
!-----------------------------------------------------------------------
               term1 = 60.51_r8
               term2 = .095_r8*wt
	       wrk   = wt*wt
               term3 = .0077_r8*wrk
               term4 = 1.61e-5_r8*wt*wrk
               term5 = (1.76_r8 + 2.52e-4_r8*wrk) * T_limitsq
               term6 = -805.89_r8 + (253.05_r8*(wt**.076_r8))
               term7 = T_limitsq
               ah    = exp( term1 - term2 + term3 - term4 - term5 + term6/term7 )
	       if( ah <= 0._r8 ) then
	          write(iulog,*) 'ratecon: ah <= 0 at i,k, = ',i,k
	          write(iulog,*) 'ratecon: term1,term2,term3,term4,term5,term6,term7,wt,T_limit,ah = ', &
	                               term1,term2,term3,term4,term5,term6,term7,wt,T_limit,ah 
	       end if

	       wrk      = .25_r8*sadsulf
               rad_sulf = max( rad_sulfate(i,k),1.e-7_r8 )
!-----------------------------------------------------------------------
!     N2O5 + H2O(liq) =>  2.00*HNO3  Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
               if( n2o5vmr > small ) then
                  term0 = -25.5265_r8 - wt*(.133188_r8 - wt*(.00930846_r8 - 9.0194e-5_r8*wt))
                  term1 = 9283.76_r8 + wt*(115.345_r8 - wt*(5.19258_r8 - .0483464_r8*wt))
                  term2 = -851801._r8 - wt*(22191.2_r8 - wt*(766.916_r8 - 6.85427_r8*wt))
                  gprob_n2o5(i,k) = exp( term0 + T_limiti*(term1 + term2*T_limiti) )
                  rxt(i,k,rid_het1) = max( 0._r8,wrk*av_n2o5*gprob_n2o5(i,k) )
               end if

!-----------------------------------------------------------------------
!     ClONO2 + H2O(liq) =  HOCl + HNO3   Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
! 	... NOTE: Aerosol radius in units of cm.
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     	... Radius sulfate set (from sad module)
!           Set min radius to 0.001 microns (1e-7 cm)
!           Typical radius is 0.1 microns (1e-5 cm)
!           f_cnt may go negative under if not set.
!-----------------------------------------------------------------------
                  C_cnt         = 1474._r8*T_limitsq
                  S_cnt         = .306_r8 + 24._r8*T_limiti
                  term1         = exp( -S_cnt*molar_h2so4 )
                  H_cnt         = 1.6e-6_r8 * exp( 4710._r8*T_limiti )*term1
                  D_cnt         = 5.e-8_r8*T_limit / vis_h2so4
                  k_h           = 1.22e12_r8*exp( -6200._r8*T_limiti )
                  k_h2o         = 1.95e10_r8*exp( -2800._r8*T_limiti )
                  k_hydr        = (k_h2o + k_h*ah)*aw
                  k_hcl         = 7.9e11_r8*ah*D_cnt*M_hcl
                  rdl_cnt       = sqrt( D_cnt/(k_hydr + k_hcl) )
                  term1         = 1._r8/tanh( rad_sulf/rdl_cnt )
                  term2         = rdl_cnt/rad_sulf
                  f_cnt         = term1 - term2
                  if( f_cnt > 0._r8 ) then
                     term1         = 4._r8*H_cnt*.082_r8*T_limit
                     term2         = sqrt( D_cnt*k_hydr )
                     Gamma_b_h2o   = term1*term2/C_cnt
                     term1         = sqrt( 1._r8 + k_hcl/k_hydr )
                     Gamma_cnt_rxn = f_cnt*Gamma_b_h2o*term1
                     Gamma_b_hcl   = Gamma_cnt_rxn*k_hcl/(k_hcl + k_hydr)
                     term1         = exp( -1374._r8*T_limiti )
                     Gamma_s       = 66.12_r8*H_cnt*M_hcl*term1
		     if( pHCl_atm > 0._r8 ) then
                        term1      = .612_r8*(Gamma_s+Gamma_b_hcl)* pCNT_atm/pHCl_atm
                        Fhcl       = 1._r8/(1._r8 + term1)
		     else
                        Fhcl       = 1._r8
		     end if
                     Gamma_s_prime     = Fhcl*Gamma_s
                     Gamma_b_hcl_prime = Fhcl*Gamma_b_hcl
                     term1         = Gamma_cnt_rxn*k_hydr
                     term2         = k_hcl + k_hydr
                     Gamma_b       = Gamma_b_hcl_prime + (term1/term2)
                     term1         = 1._r8 / (Gamma_s_prime + Gamma_b)
                     gprob_cnt     = 1._r8 / (1._r8 + term1)
                     term1         = Gamma_s_prime + Gamma_b_hcl_prime
                     term2         = Gamma_s_prime + Gamma_b
                     gprob_cnt_hcl(i,k) = gprob_cnt * term1/term2
                     gprob_cnt_h2o(i,k) = gprob_cnt - gprob_cnt_hcl(i,k)
                  else
                     gprob_cnt_h2o(i,k) = 0._r8
                     gprob_cnt_hcl(i,k) = 0._r8
                     Fhcl          = 1._r8
                  end if
                  if( clono2vmr > small ) then
                     rxt(i,k,rid_het2) = max( 0._r8,wrk*av_clono2*gprob_cnt_h2o(i,k) )
                  end if

!-----------------------------------------------------------------------
!  	... BrONO2 + H2O(liq) =  HOBr + HNO3   Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
               if( brono2vmr > small ) then
                  h1    = 29.24_r8
                  h2    = -.396_r8
                  h3    = .114_r8
                  alpha = .805_r8
                  gprob_rxn = exp( h1 + h2*wt ) + h3
                  term1     = 1._r8/alpha
                  term2     = 1._r8/gprob_rxn
                  gprob_bnt_h2o(i,k) = 1._r8 / (term1 + term2)
                  rxt(i,k,rid_het3) = max( 0._r8,wrk*av_brono2*gprob_bnt_h2o(i,k) )
               end if

!-----------------------------------------------------------------------
!     	... ClONO2 + HCl(liq) =  Cl2  + HNO3  Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
               if( hclvmr > small .and. clono2vmr > small ) then
                 if ( hclvmr > clono2vmr ) then
                    rxt(i,k,rid_het4) = max( 0._r8,wrk*av_clono2*gprob_cnt_hcl(i,k) )*hcldeni
                 else
                    rxt(i,k,rid_het4) = max( 0._r8,wrk*av_clono2*gprob_cnt_hcl(i,k) )*cntdeni
                 end if
               end if

!-----------------------------------------------------------------------
!     	... HOCl + HCl(liq) =  Cl2 + H2O   Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
               if( hclvmr > small .and. hoclvmr > small ) then
!-----------------------------------------------------------------------
!     	... Radius sulfate set (from sad module)
!           Set min radius to 0.001 microns (1e-7 cm)
!           Typical radius is 0.1 microns (1e-5 cm)
!           f_hocl may go negative under if not set.
!-----------------------------------------------------------------------
	          if( pCNT_atm > 0._r8 ) then
                     D_hocl          = 6.4e-8_r8*T_limit/vis_h2so4
                     k_hocl_hcl      = 1.25e9_r8*ah*D_hocl*M_hcl
                     C_hocl          = 2009._r8*T_limitsq
                     S_hocl          = .0776_r8 + 59.18_r8*T_limiti
                     term1           = exp( -S_hocl*molar_h2so4 )
                     H_hocl          = 1.91e-6_r8 * exp( 5862.4_r8*T_limiti )*term1
                     term1           = 4._r8*H_hocl*.082_r8*T_limit
                     term2           = sqrt( D_hocl*k_hocl_hcl )
                     Gamma_hocl_rxn  = term1*term2/C_hocl
                     rdl_hocl        = sqrt( D_hocl/k_hocl_hcl )
                     term1           = 1._r8/tanh( rad_sulf/rdl_hocl )
                     term2           = rdl_hocl/rad_sulf
                     f_hocl          = term1 - term2
                     if( f_hocl > 0._r8 ) then
                        term1           = 1._r8 / (f_hocl*Gamma_hocl_rxn*Fhcl)
                        gprob_hocl_hcl(i,k)  = 1._r8 / (1._r8 + term1)
                     else
                        gprob_hocl_hcl(i,k)  = 0._r8
                     end if

                     if ( hclvmr > hoclvmr ) then
                       rxt(i,k,rid_het5) = max( 0._r8,wrk*av_hocl*gprob_hocl_hcl(i,k) )*hcldeni
                     else
                       rxt(i,k,rid_het5) = max( 0._r8,wrk*av_hocl*gprob_hocl_hcl(i,k) )*hocldeni
                     end if

	          end if
               end if

!-----------------------------------------------------------------------
!     	... HOBr + HCl(liq) =  BrCl + H2O  Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
               if( hclvmr > small .and. hobrvmr > small ) then
!-----------------------------------------------------------------------
!   	... Radius sulfate set (from sad module)
!           Set min radius to 0.001 microns (1e-7 cm)
!           Typical radius is 0.1 microns (1e-5 cm)
!           f_hobr may go negative under if not set.
!-----------------------------------------------------------------------
                  C_hobr          = 1477._r8*T_limitsq
                  D_hobr          = 9.e-9_r8
!-----------------------------------------------------------------------
!     	...  Taken from Waschewsky and Abbat
!            Dave Hanson (PC) suggested we divide this rc by eight to agee
!            with his data (Hanson, in press, 2002).
!            k1=k2*Mhcl for gamma(HOBr)
!-----------------------------------------------------------------------
                  k_wasch         = .125_r8 * exp( .542_r8*wt - 6440._r8*T_limiti + 10.3_r8)
!-----------------------------------------------------------------------
!     	... Taken from Hanson 2002.
!-----------------------------------------------------------------------
                  H_hobr          = exp( -9.86_r8 + 5427._r8*T_limiti )
                  k_dl            = 7.5e14_r8*D_hobr*2._r8                        ! or  7.5e14*D *(2nm)
!-----------------------------------------------------------------------
!  	... If k_wasch is GE than the diffusion limit...
!-----------------------------------------------------------------------
		  if( M_hcl > 0._r8 ) then
                     if( k_wasch >= k_dl ) then
                        k_hobr_hcl   = k_dl * M_hcl
                     else
                        k_hobr_hcl   = k_wasch * M_hcl
                     end if
                     term1           = 4._r8*H_hobr*.082_r8*T_limit
                     term2           = sqrt( D_hobr*k_hobr_hcl )
                     tmp             = rad_sulf/term2
                     Gamma_hobr_rxn  = term1*term2/C_hobr
                     rdl_hobr        = sqrt( D_hobr/k_hobr_hcl )
		     if( tmp < 1.e2_r8 ) then
                        term1           = 1._r8/tanh( rad_sulf/rdl_hobr )
		     else
                        term1           = 1._r8
		     end if
                     term2           = rdl_hobr/rad_sulf
                     f_hobr          = term1 - term2
                     if( f_hobr > 0._r8 ) then
                        term1            = 1._r8 / (f_hobr*Gamma_hobr_rxn)
                        gprob_hobr_hcl(i,k)   = 1._r8 / (1._r8 + term1)
                     else
                         gprob_hobr_hcl(i,k)  = 0._r8
                     end if

                     if ( hclvmr > hobrvmr ) then
                        rxt(i,k,rid_het6) = max( 0._r8,wrk*av_hobr*gprob_hobr_hcl(i,k) )*hcldeni
                     else
                        rxt(i,k,rid_het6) = max( 0._r8,wrk*av_hobr*gprob_hobr_hcl(i,k) )*hobrdeni    
                     end if           

		  end if
               end if
            end if has_sadsulf

has_sadnat : &
	    if( sadnat > 0._r8 ) then
	       wrk = .25_r8*sadnat
!-----------------------------------------------------------------------
!     	... N2O5 + H2O(s) => 2HNO3  NAT Aerosol Reaction
!-----------------------------------------------------------------------
               if( n2o5vmr > small ) then
!-----------------------------------------------------------------------
!     ... gprob based on JPL06-2 for NAT.
!         also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
!                 gprob_tot     = 4.e-4
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het7)  = wrk*av_n2o5*4.e-4_r8
	       end if

!-----------------------------------------------------------------------
!     ClONO2 + H2O(s) => HNO3 + HOCl  NAT Aerosol Reaction
!-----------------------------------------------------------------------
               if( clono2vmr > small ) then
!-----------------------------------------------------------------------
!     ... gprob based on JPL06-2 for NAT.
!         also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
!                 gprob_tot    = 0.004
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het8) = wrk*av_clono2*4.0e-3_r8
   	       end if

!-----------------------------------------------------------------------
!     	... ClONO2 + HCl(s) => HNO3 + Cl2, NAT Aerosol Reaction
!-----------------------------------------------------------------------
               if( hclvmr > small ) then
                  if( clono2vmr > small ) then
!-----------------------------------------------------------------------
!     ... gprob based on JPL06-2 for NAT.
!         also see Hanson and Ravi, JPC, 96, 2682-2691, 1992.
!                 gprob_tot   = 0.2
!-----------------------------------------------------------------------
                     if ( hclvmr > clono2vmr ) then
                        rxt(i,k,rid_het9) = wrk*av_clono2*0.2_r8*hcldeni
                     else
                        rxt(i,k,rid_het9) = wrk*av_clono2*0.2_r8*cntdeni  
                     end if
                  end if

!-----------------------------------------------------------------------
!     	... HOCl + HCl(s) => H2O + Cl2  NAT Aerosol Reaction
!-----------------------------------------------------------------------
                  if( hoclvmr > small ) then
!-----------------------------------------------------------------------
!     ... gprob based on JPL06-2 for NAT.
!         see Hanson and Ravi, JPC, 96, 2682-2691, 1992.
!         and      Abbatt and Molina, GRL, 19, 461-464, 1992.
!                 gprob_tot   = 0.1
!-----------------------------------------------------------------------
                     if ( hclvmr > hoclvmr ) then
                        rxt(i,k,rid_het10) = wrk*av_hocl*0.1_r8*hcldeni
                     else
                        rxt(i,k,rid_het10) = wrk*av_hocl*0.1_r8*hocldeni
                     end if
                  end if
               end if

!-----------------------------------------------------------------------
!     	... BrONO2 + H2O(s) => HOBr + HNO3  NAT Aerosol Reaction
!-----------------------------------------------------------------------
               if( brono2vmr > small ) then
!-----------------------------------------------------------------------
!       ... Personel Communication, 11/4/99, David Hanson
!                 gprob_tot   = 0.3
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het11) = wrk*av_brono2*0.3_r8
               end if
            end if has_sadnat

has_sadice : &
	    if( sadice > 0._r8 ) then
	       wrk = .25_r8*sadice
!-----------------------------------------------------------------------
!     N2O5 + H2O(s) => 2HNO3  ICE Aerosol Reaction
!-----------------------------------------------------------------------
               if( n2o5vmr > small ) then
!-----------------------------------------------------------------------
!       ... gprob based on JPL06-2 for ICE.
!           also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
!                 gprob_tot    = .02
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het12) = wrk*av_n2o5*0.02_r8
 	       end if
!-----------------------------------------------------------------------
!     	... ClONO2 + H2O(s) => HNO3 + HOCl  ICE Aerosol Reaction
!-----------------------------------------------------------------------
               if( clono2vmr > small ) then
!-----------------------------------------------------------------------
!     	... gprob based on JPL06-2 for ICE.
!     	    also see Hanson and Ravi, JGR, 96, 17307-17314, 1991.
!                 gprob_tot    = .3
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het13) = wrk*av_clono2*0.3_r8
	       end if

!-----------------------------------------------------------------------
!     	... BrONO2 + H2O(s) => HNO3 + HOBr  ICE Aerosol Reaction
!-----------------------------------------------------------------------
               if( brono2vmr > small ) then
!-----------------------------------------------------------------------
!     	... gprob based on JPL06-2 for ICE.
!           also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
!           could be as high as 1.0
!                 gprob_tot    = .3
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het14) = wrk*av_brono2*0.3_r8
      	       end if

!-----------------------------------------------------------------------
!     ClONO2 + HCl(s) => HNO3 + Cl2, ICE Aerosol Reaction
!-----------------------------------------------------------------------
               if( hclvmr > small ) then
                  if( clono2vmr > small ) then
!-----------------------------------------------------------------------
!       ... gprob based on JPL06-2 for ICE.
!           also see Hanson and Ravi, GRL, 15, 17-20, 1988.
!           also see Lue et al.,
!                 gprob_tot    = .3
!-----------------------------------------------------------------------
                     if ( hclvmr > clono2vmr ) then
                        rxt(i,k,rid_het15) = wrk*av_clono2*0.3_r8*hcldeni
                     else
                        rxt(i,k,rid_het15) = wrk*av_clono2*0.3_r8*cntdeni
                     end if

                  end if
!
!-----------------------------------------------------------------------
!     	... HOCl + HCl(s) => H2O + Cl2, ICE Aerosol Reaction
!-----------------------------------------------------------------------
                  if( hoclvmr > small .and. hclvmr > small ) then
!-----------------------------------------------------------------------
!       ... gprob based on JPL06-2 for ICE.
!           also see Hanson and Ravi, JPC, 96, 2682-2691, 1992.
!           also see Abbatt and Molina, GRL, 19, 461-464, 1992.
!                 gprob_tot   = .2
!-----------------------------------------------------------------------
                     if ( hclvmr > hoclvmr ) then
                        rxt(i,k,rid_het16) = wrk*av_hocl*0.2_r8*hcldeni
                     else
                        rxt(i,k,rid_het16) = wrk*av_hocl*0.2_r8*hocldeni
                     end if

                  end if

!-----------------------------------------------------------------------
!     HOBr + HCl(s) => H2O + BrCl, ICE Aerosol Reaction
!-----------------------------------------------------------------------
                  if( hobrvmr > small .and. hclvmr > small ) then
!-----------------------------------------------------------------------
!       ... gprob based on JPL06-2 for ICE.
!           Abbatt GRL, 21, 665-668, 1994.
!                    gprob_tot   = .3
!-----------------------------------------------------------------------
                    if ( hclvmr > hobrvmr ) then
                       rxt(i,k,rid_het17) = wrk*av_hobr*0.3_r8*hcldeni
                    else
                       rxt(i,k,rid_het17) = wrk*av_hobr*0.3_r8*hobrdeni
                    end if
                  end if
	       end if
            end if has_sadice
         end do column_loop
      end do Level_loop

      end subroutine ratecon_sfstrat

      end module mo_strato_rates
